Data analysis for:

Why are telomeres the length that they are? Insight from a phylogenetic comparative analysis

Workflow of the analyses produced by Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.


Library

library(AICcmodavg)
library(ape)
library(caper)
library(car)
library(coda)
library(extrafont)
library(geiger)
library(graph)
library(kableExtra)
library(MuMIn)
library(nlme)
library(pbapply)
library(phylopath)
library(phytools)
library(plotrix)
library(rphylopic)
library(scales)
library(sensiPhy)
library(shape)
library(xtable)
R.version
>                _                           
> platform       aarch64-apple-darwin20      
> arch           aarch64                     
> os             darwin20                    
> system         aarch64, darwin20           
> status                                     
> major          4                           
> minor          4.1                         
> year           2024                        
> month          06                          
> day            14                          
> svn rev        86737                       
> language       R                           
> version.string R version 4.4.1 (2024-06-14)
> nickname       Race for Your Life
sessionInfo()
> R version 4.4.1 (2024-06-14)
> Platform: aarch64-apple-darwin20
> Running under: macOS Sonoma 14.6.1
> 
> Matrix products: default
> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
> 
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
> 
> time zone: America/Phoenix
> tzcode source: internal
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base     
> 
> other attached packages:
>  [1] rmarkdown_2.27       xtable_1.8-4         shape_1.4.6.1        sensiPhy_0.8.5       ggplot2_3.5.1        phylolm_2.6.2       
>  [7] scales_1.3.0         rphylopic_1.4.0.9000 plotrix_3.8-4        phylopath_1.3.0      pbapply_1.7-2        nlme_3.1-165        
> [13] MuMIn_1.48.4         kableExtra_1.4.0     graph_1.82.0         BiocGenerics_0.50.0  geiger_2.0.11        phytools_2.3-0      
> [19] maps_3.4.2           extrafont_0.19       coda_0.19-4.1        car_3.1-2            carData_3.0-5        caper_1.0.3         
> [25] mvtnorm_1.2-5        MASS_7.3-60.2        ape_5.8              AICcmodavg_2.3-3    
> 
> loaded via a namespace (and not attached):
>   [1] mnormt_2.1.1            gridExtra_2.3           phangorn_2.11.1         rlang_1.1.4             magrittr_2.0.3         
>   [6] compiler_4.4.1          png_0.1-8               systemfonts_1.1.0       vctrs_0.6.5             combinat_0.0-8         
>  [11] quadprog_1.5-8          stringr_1.5.1           pkgconfig_2.0.3         fastmap_1.2.0           ggraph_2.2.1           
>  [16] labeling_0.4.3          utf8_1.2.4              subplex_1.9             deSolve_1.40            purrr_1.0.2            
>  [21] xfun_0.45               cachem_1.1.0            clusterGeneration_1.3.8 jsonlite_1.8.8          highr_0.11             
>  [26] tweenr_2.0.3            jpeg_0.1-10             VGAM_1.1-11             parallel_4.4.1          R6_2.5.1               
>  [31] bslib_0.7.0             stringi_1.8.4           parallelly_1.38.0       extrafontdb_1.0         jquerylib_0.1.4        
>  [36] numDeriv_2016.8-1.1     Rcpp_1.0.12             iterators_1.0.14        knitr_1.47              future.apply_1.11.2    
>  [41] optimParallel_1.0-2     ggm_2.5.1               base64enc_0.1-3         Matrix_1.7-0            splines_4.4.1          
>  [46] igraph_2.0.3            tidyselect_1.2.1        viridis_0.6.5           yaml_2.3.8              rstudioapi_0.16.0      
>  [51] abind_1.4-5             doParallel_1.0.17       codetools_0.2-20        curl_5.2.1              listenv_0.9.1          
>  [56] lattice_0.22-6          tibble_3.2.1            withr_3.0.0             evaluate_0.24.0         future_1.34.0          
>  [61] survival_3.6-4          polyclip_1.10-7         xml2_1.3.6              BiocManager_1.30.23     pillar_1.9.0           
>  [66] foreach_1.5.2           stats4_4.4.1            generics_0.1.3          grImport2_0.3-1         munsell_0.5.1          
>  [71] globals_0.16.3          glue_1.7.0              unmarked_1.4.1          scatterplot3d_0.3-44    tools_4.4.1            
>  [76] graphlayouts_1.1.1      XML_3.99-0.17           tidygraph_1.3.1         fastmatch_1.1-4         grid_4.4.1             
>  [81] tidyr_1.3.1             Rttf2pt1_1.3.12         colorspace_2.1-0        ggforce_0.4.2           cli_3.6.3              
>  [86] DEoptim_2.2-8           rsvg_2.6.0              fansi_1.0.6             expm_0.999-9            viridisLite_0.4.2      
>  [91] svglite_2.1.3           dplyr_1.1.4             gtable_0.3.5            sass_0.4.9              digest_0.6.36          
>  [96] ggrepel_0.9.5           farver_2.1.2            memoise_2.0.1           htmltools_0.5.8.1       lifecycle_1.0.4        
> [101] httr_1.4.7
set.seed(80)

## Dataset

data <- read.csv('../TRF-anaysis-no-domesticated.csv')
str(data)
> 'data.frame': 131 obs. of  27 variables:
>  $ Common.name               : chr  "European seabass" "Skipjack trevally" "Green swordtail" "Southern platyfish" ...
>  $ Scientific_name           : chr  "Dicentrarchus_labrax" "Pseudocaranx_wrighti" "Xiphophorus_hellerii" "Xiphophorus_maculatus" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Fish" "Fish" "Fish" "Fish" ...
>  $ Order                     : chr  "Moroniformes" "Carangiformes" "Cyprinodontiformes" "Cyprinodontiformes" ...
>  $ Family                    : chr  "Moronidae" "Carangidae" "Poeciliidae" "Poeciliidae" ...
>  $ Genus                     : chr  "Dicentrarchus" "Pseudocaranx" "Xiphophorus" "Xiphophorus" ...
>  $ Species                   : chr  "labrax" "wrighti" "hellerii" "maculatus" ...
>  $ Endo_ectotherm            : chr  "ecto" "ecto" "ecto" "ecto" ...
>  $ Adult_mass_grams          : num  6600 3000 3.7 1.8 1600 ...
>  $ log_mass                  : num  3.82 3.477 0.568 0.255 3.204 ...
>  $ Lifespan_years            : num  15 13 5 5 41 47 18 10 7 20 ...
>  $ Loglifespan               : num  1.176 1.114 0.699 0.699 1.613 ...
>  $ Average_Telomere_Length_kb: num  3.44 4.06 4.25 4.25 4.6 4.62 5.26 5.48 5.56 6.14 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
>  $ Tissue.type.for.TA        : chr  NA NA NA NA ...
>  $ Source.for.TA             : chr  NA NA NA NA ...
>  $ Tissue.type.for.TL        : chr  "Blood" "Muscle" "Multiple somatic tissues" "Multiple somatic tissues" ...
>  $ Methodology.for.TL        : chr  "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ TRF_type                  : chr  "A" "A" "A" "A" ...
>  $ Source.for.TL             : chr  "Horn et al. 2008" "Izzo 2010" "Downs et al. 2012" "Downs et al. 2012" ...
>  $ Source.for.Lifespan       : chr  "AnAge" "Ratnasingham and Herbert 2007" "AnAge" "National park aquarium" ...
>  $ Source.for.Mass           : chr  "AnAge" "Marinewise" "Meyer et al. 2006" "McKenzie et al. 1983" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
head(data)
>          Common.name       Scientific_name    Domain Kingdom   Phylum Class              Order      Family         Genus   Species
> 1   European seabass  Dicentrarchus_labrax Eukaryota Metazoa Chordata  Fish       Moroniformes   Moronidae Dicentrarchus    labrax
> 2  Skipjack trevally  Pseudocaranx_wrighti Eukaryota Metazoa Chordata  Fish      Carangiformes  Carangidae  Pseudocaranx   wrighti
> 3    Green swordtail  Xiphophorus_hellerii Eukaryota Metazoa Chordata  Fish Cyprinodontiformes Poeciliidae   Xiphophorus  hellerii
> 4 Southern platyfish Xiphophorus_maculatus Eukaryota Metazoa Chordata  Fish Cyprinodontiformes Poeciliidae   Xiphophorus maculatus
> 5           Goldfish     Carassius_auratus Eukaryota Metazoa Chordata  Fish      Cypriniformes  Cyprinidae     Carassius   auratus
> 6        Common carp       Cyprinus_carpio Eukaryota Metazoa Chordata  Fish      Cypriniformes  Cyprinidae      Cyprinus    carpio
>   Endo_ectotherm Adult_mass_grams  log_mass Lifespan_years Loglifespan Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA
> 1           ecto           6600.0 3.8195439             15    1.176091                       3.44                  NA               <NA>
> 2           ecto           3000.0 3.4771213             13    1.113943                       4.06                  NA               <NA>
> 3           ecto              3.7 0.5682017              5    0.698970                       4.25                  NA               <NA>
> 4           ecto              1.8 0.2552725              5    0.698970                       4.25                  NA               <NA>
> 5           ecto           1600.0 3.2041200             41    1.612784                       4.60                  NA               <NA>
> 6           ecto          10600.0 4.0253059             47    1.672098                       4.62                  NA               <NA>
>   Source.for.TA       Tissue.type.for.TL Methodology.for.TL TRF_type     Source.for.TL           Source.for.Lifespan      Source.for.Mass
> 1          <NA>                    Blood    TRF (denatured)        A  Horn et al. 2008                         AnAge                AnAge
> 2          <NA>                   Muscle    TRF (denatured)        A         Izzo 2010 Ratnasingham and Herbert 2007           Marinewise
> 3          <NA> Multiple somatic tissues    TRF (denatured)        A Downs et al. 2012                         AnAge    Meyer et al. 2006
> 4          <NA> Multiple somatic tissues    TRF (denatured)        A Downs et al. 2012        National park aquarium McKenzie et al. 1983
> 5          <NA>                   Muscle    TRF (denatured)        A         Izzo 2010                         AnAge             Fishbase
> 6          <NA>                   Muscle    TRF (denatured)        A         Izzo 2010                         AnAge                AnAge
>   Captivity Notes
> 1         0      
> 2         0      
> 3         0      
> 4         0      
> 5         0      
> 6         0
unique(data$Class)
> [1] "Fish"     "Reptilia" "Mammalia" "Aves"
dat <- data
names(dat)
>  [1] "Common.name"                "Scientific_name"            "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                      "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                    "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"             "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "Telomerase_activity"        "Tissue.type.for.TA"         "Source.for.TA"              "Tissue.type.for.TL"        
> [21] "Methodology.for.TL"         "TRF_type"                   "Source.for.TL"              "Source.for.Lifespan"       
> [25] "Source.for.Mass"            "Captivity"                  "Notes"
unique(is.na(dat$log_mass))
> [1] FALSE
## Trees

##write.table(data$Scientific_name, 'spp.tree.08-03-24.txt', row.names = FALSE, col.names = FALSE)

full_data_tree <- read.tree("spp.tree.08-03-24.nwk")
is.ultrametric(full_data_tree)
> [1] FALSE
full_data_tree <- force.ultrametric(full_data_tree)
> ***************************************************************
> *                          Note:                              *
> *    force.ultrametric does not include a formal method to    *
> *    ultrametricize a tree & should only be used to coerce    *
> *   a phylogeny that fails is.ultrametric due to rounding --  *
> *    not as a substitute for formal rate-smoothing methods.   *
> ***************************************************************
is.ultrametric(full_data_tree)
> [1] TRUE
full_data_tree
> 
> Phylogenetic tree with 138 tips and 137 internal nodes.
> 
> Tip labels:
>   Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
>   , '14', '12', '48', '57', '49', ...
> 
> Rooted; includes branch lengths.
rownames(dat) <- dat$Scientific_name


## Full_tree

check <- name.check(full_data_tree, dat)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pruned_data_tree <- drop.tip(full_data_tree, rm_phy)
pruned_dat <- subset(dat, subset = dat$Scientific_name %in% pruned_data_tree$tip, select = names(dat))
str(pruned_dat)
> 'data.frame': 123 obs. of  27 variables:
>  $ Common.name               : chr  "European seabass" "Skipjack trevally" "Green swordtail" "Southern platyfish" ...
>  $ Scientific_name           : chr  "Dicentrarchus_labrax" "Pseudocaranx_wrighti" "Xiphophorus_hellerii" "Xiphophorus_maculatus" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Fish" "Fish" "Fish" "Fish" ...
>  $ Order                     : chr  "Moroniformes" "Carangiformes" "Cyprinodontiformes" "Cyprinodontiformes" ...
>  $ Family                    : chr  "Moronidae" "Carangidae" "Poeciliidae" "Poeciliidae" ...
>  $ Genus                     : chr  "Dicentrarchus" "Pseudocaranx" "Xiphophorus" "Xiphophorus" ...
>  $ Species                   : chr  "labrax" "wrighti" "hellerii" "maculatus" ...
>  $ Endo_ectotherm            : chr  "ecto" "ecto" "ecto" "ecto" ...
>  $ Adult_mass_grams          : num  6600 3000 3.7 1.8 1600 ...
>  $ log_mass                  : num  3.82 3.477 0.568 0.255 3.204 ...
>  $ Lifespan_years            : num  15 13 5 5 41 47 18 10 7 20 ...
>  $ Loglifespan               : num  1.176 1.114 0.699 0.699 1.613 ...
>  $ Average_Telomere_Length_kb: num  3.44 4.06 4.25 4.25 4.6 4.62 5.26 5.48 5.56 6.14 ...
>  $ Telomerase_activity       : int  NA NA NA NA NA NA NA NA NA NA ...
>  $ Tissue.type.for.TA        : chr  NA NA NA NA ...
>  $ Source.for.TA             : chr  NA NA NA NA ...
>  $ Tissue.type.for.TL        : chr  "Blood" "Muscle" "Multiple somatic tissues" "Multiple somatic tissues" ...
>  $ Methodology.for.TL        : chr  "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ TRF_type                  : chr  "A" "A" "A" "A" ...
>  $ Source.for.TL             : chr  "Horn et al. 2008" "Izzo 2010" "Downs et al. 2012" "Downs et al. 2012" ...
>  $ Source.for.Lifespan       : chr  "AnAge" "Ratnasingham and Herbert 2007" "AnAge" "National park aquarium" ...
>  $ Source.for.Mass           : chr  "AnAge" "Marinewise" "Meyer et al. 2006" "McKenzie et al. 1983" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
head(pruned_dat)
>                              Common.name       Scientific_name    Domain Kingdom   Phylum Class              Order      Family
> Dicentrarchus_labrax    European seabass  Dicentrarchus_labrax Eukaryota Metazoa Chordata  Fish       Moroniformes   Moronidae
> Pseudocaranx_wrighti   Skipjack trevally  Pseudocaranx_wrighti Eukaryota Metazoa Chordata  Fish      Carangiformes  Carangidae
> Xiphophorus_hellerii     Green swordtail  Xiphophorus_hellerii Eukaryota Metazoa Chordata  Fish Cyprinodontiformes Poeciliidae
> Xiphophorus_maculatus Southern platyfish Xiphophorus_maculatus Eukaryota Metazoa Chordata  Fish Cyprinodontiformes Poeciliidae
> Carassius_auratus               Goldfish     Carassius_auratus Eukaryota Metazoa Chordata  Fish      Cypriniformes  Cyprinidae
> Cyprinus_carpio              Common carp       Cyprinus_carpio Eukaryota Metazoa Chordata  Fish      Cypriniformes  Cyprinidae
>                               Genus   Species Endo_ectotherm Adult_mass_grams  log_mass Lifespan_years Loglifespan
> Dicentrarchus_labrax  Dicentrarchus    labrax           ecto           6600.0 3.8195439             15    1.176091
> Pseudocaranx_wrighti   Pseudocaranx   wrighti           ecto           3000.0 3.4771213             13    1.113943
> Xiphophorus_hellerii    Xiphophorus  hellerii           ecto              3.7 0.5682017              5    0.698970
> Xiphophorus_maculatus   Xiphophorus maculatus           ecto              1.8 0.2552725              5    0.698970
> Carassius_auratus         Carassius   auratus           ecto           1600.0 3.2041200             41    1.612784
> Cyprinus_carpio            Cyprinus    carpio           ecto          10600.0 4.0253059             47    1.672098
>                       Average_Telomere_Length_kb Telomerase_activity Tissue.type.for.TA Source.for.TA       Tissue.type.for.TL
> Dicentrarchus_labrax                        3.44                  NA               <NA>          <NA>                    Blood
> Pseudocaranx_wrighti                        4.06                  NA               <NA>          <NA>                   Muscle
> Xiphophorus_hellerii                        4.25                  NA               <NA>          <NA> Multiple somatic tissues
> Xiphophorus_maculatus                       4.25                  NA               <NA>          <NA> Multiple somatic tissues
> Carassius_auratus                           4.60                  NA               <NA>          <NA>                   Muscle
> Cyprinus_carpio                             4.62                  NA               <NA>          <NA>                   Muscle
>                       Methodology.for.TL TRF_type     Source.for.TL           Source.for.Lifespan      Source.for.Mass Captivity Notes
> Dicentrarchus_labrax     TRF (denatured)        A  Horn et al. 2008                         AnAge                AnAge         0      
> Pseudocaranx_wrighti     TRF (denatured)        A         Izzo 2010 Ratnasingham and Herbert 2007           Marinewise         0      
> Xiphophorus_hellerii     TRF (denatured)        A Downs et al. 2012                         AnAge    Meyer et al. 2006         0      
> Xiphophorus_maculatus    TRF (denatured)        A Downs et al. 2012        National park aquarium McKenzie et al. 1983         0      
> Carassius_auratus        TRF (denatured)        A         Izzo 2010                         AnAge             Fishbase         0      
> Cyprinus_carpio          TRF (denatured)        A         Izzo 2010                         AnAge                AnAge         0
pruned_data_tree
> 
> Phylogenetic tree with 123 tips and 122 internal nodes.
> 
> Tip labels:
>   Urolophus_paucimaculatus, Callorhinchus_milii, Oncorhynchus_mykiss, Lutjanus_argentimaculatus, Acanthopagrus_schlegelii, Dicentrarchus_labrax, ...
> Node labels:
>   , '14', '12', '48', '57', '49', ...
> 
> Rooted; includes branch lengths.
name.check(pruned_data_tree, pruned_dat)
> [1] "OK"
## Checking the distribution of lifespan

plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()

par(new = TRUE)
hist(pruned_dat$Lifespan_years, main = "Raw variable", las = 1, xlab = 'Lifespan (years)')
box()

plot(NA, xlim = c(0, 200), ylim = c(0, 60), type = 'n', las = 1, xlab = '', ylab = '', axes = FALSE)
grid()

par(new = TRUE)
hist(log1p(pruned_dat$Lifespan_years), main = "Log-transformed", xlab = 'Log lifespan (years)', las = 1)
box()

## Log-transforming for a better fit

pruned_dat$log.lifespan <- log1p(pruned_dat$Lifespan_years)
pruned_dat$log.mass <- log1p(pruned_dat$Adult_mass_grams)



Implementing phylogenetic path analysis as suggested by the editor


## Defining models

path.mod <- pruned_dat
names(path.mod)
>  [1] "Common.name"                "Scientific_name"            "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                      "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                    "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"             "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "Telomerase_activity"        "Tissue.type.for.TA"         "Source.for.TA"              "Tissue.type.for.TL"        
> [21] "Methodology.for.TL"         "TRF_type"                   "Source.for.TL"              "Source.for.Lifespan"       
> [25] "Source.for.Mass"            "Captivity"                  "Notes"                      "log.lifespan"              
> [29] "log.mass"
path.mod$Average_Telomere_Length_kb <- log(path.mod$Average_Telomere_Length_kb)

names(path.mod)[c(16, 11, 29, 28)] <- c('TL', 'TB', 'BM', 'LS')

models <- define_model_set(
    "null" = c(TL ~ TL, TB ~ LS, TB ~ BM, LS ~ BM),
    "(a)" = c(TL ~ LS, TB ~ LS, TB ~ BM, BM ~ LS),
    "(b)" = c(TL ~ BM, TB ~ LS, TB ~ BM, BM ~ LS),
    "(c)" = c(TL ~ TB, TB ~ LS, TB ~ BM, BM ~ LS),
    "(d)" = c(TL ~ BM + LS + TB),
    "(e)" = c(TL ~ LS),
    "(f)" = c(TL ~ BM),
    "(g)" = c(TL ~ TB))


#png("figureS1-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS1-revised.pdf")

plot_model_set(models, nrow = 4) + ggplot2::scale_y_continuous(expand = c(0.1, 0)) + theme(plot.margin=unit(c(0.5, 3, 0.5, 3), "cm"))

#dev.off()

mo <- phylo_path(models, path.mod, pruned_data_tree)
s <- summary(mo)
theme_set(theme_bw())

#png("figureS2-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS2-revised.pdf")

bar.pl <- tapply(s$w, list(s$model), print)
> [1] 0.06045404
> [1] 0.03543754
> [1] 0.9037687
> [1] 3.83462e-08
> [1] 3.006297e-11
> [1] 9.501932e-12
> [1] 9.928878e-12
> [1] 0.0003397325
p.val <- s[names(bar.pl), ]

barplot(NA, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), axes = FALSE)

grid()

par(new = TRUE)

barplot(bar.pl, beside = TRUE, space = FALSE, horiz = TRUE, xlim = c(0, 1), yaxt = 'n', ylab = 'Candidate models', xlab = 'Weight of evidence (w)', col = c(rep('gray', 2), alpha('red', 0.3), rep('gray', 5)))

axis(2, at = (0:7)+0.5, labels = names(bar.pl), las = 1)
box()

text(x = round(bar.pl, 3)+0.05, y = (0:7)+0.5, labels = round(p.val$p, 3))
legend('topright', legend = 'within 2 CICc', fil = alpha('red', 0.3), bty = 'n')

#dev.off()

pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model0)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   235.1572 252.0303 -111.5786
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.6254488 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.629715 0.3916907  9.266787  0.0000
> Endo_ectothermendo -0.819368 0.2207192 -3.712263  0.0003
> log.lifespan       -0.134538 0.1130636 -1.189932  0.2364
> log.mass           -0.033764 0.0258447 -1.306413  0.1939
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.041              
> log.lifespan       -0.502 -0.101       
> log.mass            0.075  0.013 -0.678
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -2.34917752 -0.02875553  0.50192730  0.92152006  2.39155282 
> 
> Residual standard error: 0.8117616 
> Degrees of freedom: 123 total; 119 residual
brown.model <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corBrownian(phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(brown.model)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   272.3855 286.4464 -131.1927
> 
> Correlation Structure: corBrownian
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> numeric(0)
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.901248 0.8933978  4.366753  0.0000
> Endo_ectothermendo -1.172484 0.2869553 -4.085949  0.0001
> log.lifespan       -0.142480 0.1268266 -1.123421  0.2635
> log.mass           -0.048054 0.0328901 -1.461046  0.1466
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.025              
> log.lifespan       -0.268 -0.136       
> log.mass           -0.086  0.108 -0.499
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.11576421  0.04130495  0.29665511  0.48025508  1.15040877 
> 
> Residual standard error: 1.926538 
> Degrees of freedom: 123 total; 119 residual
VerTree <- pruned_data_tree
VerTree$edge.length <- pruned_data_tree$edge.length * 100

OU.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corMartins(1, phy = VerTree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(OU.model1)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   283.5872 300.4603 -135.7936
> 
> Correlation Structure: corMartins
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> alpha 
>     1 
> 
> Coefficients:
>                         Value  Std.Error   t-value p-value
> (Intercept)         3.0127181 0.27461330 10.970765  0.0000
> Endo_ectothermendo -0.0099686 0.15910482 -0.062654  0.9501
> log.lifespan       -0.1146562 0.11583847 -0.989794  0.3243
> log.mass           -0.0007418 0.02556906 -0.029011  0.9769
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.108              
> log.lifespan       -0.748 -0.291       
> log.mass            0.131  0.098 -0.650
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.2551171 -0.5625105 -0.0212407  0.5130302  1.9228926 
> 
> Residual standard error: 0.7298431 
> Degrees of freedom: 123 total; 119 residual
## Model selection procedure based on AIC criterion

output.models <- model.sel(brown.model, pagel.model0, OU.model1)

sel.table <- round(as.data.frame(output.models)[6:10], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table
>                     Model K   logLik    AICc  delta weight
> pagel.model0 pagel.model0 6 -111.579 235.881  0.000      1
> brown.model   brown.model 5 -131.193 272.898 37.017      0
> OU.model1       OU.model1 6 -135.794 284.311 48.430      0
anova(pagel.model0)
> Denom. DF: 119 
>                numDF  F-value p-value
> (Intercept)        1 84.64587  <.0001
> Endo_ectotherm     1 16.64769  0.0001
> log.lifespan       1  7.96289  0.0056
> log.mass           1  1.70672  0.1939
pagel.mod.red <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.mod.red)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   234.9012 248.9622 -112.4506
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.6138179 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.662657 0.3875853  9.449939  0.0000
> Endo_ectothermendo -0.809681 0.2203418 -3.674661  0.0004
> log.lifespan       -0.233895 0.0830949 -2.814794  0.0057
> 
>  Correlation: 
>                    (Intr) End_ct
> Endo_ectothermendo -0.043       
> log.lifespan       -0.620 -0.126
> 
> Standardized residuals:
>          Min           Q1          Med           Q3          Max 
> -2.293252813 -0.003891462  0.546445007  0.947511030  2.262724689 
> 
> Residual standard error: 0.8098694 
> Degrees of freedom: 123 total; 120 residual
anova(pagel.model0, pagel.mod.red) # no difference between the models
>               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> pagel.model0      1  6 235.1572 252.0303 -111.5786                        
> pagel.mod.red     2  5 234.9013 248.9622 -112.4506 1 vs 2 1.744052  0.1866
mod.plot <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm*log.lifespan, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML") ## for plotting purposes. Notice the slopes are not significant, but the main effects are!
summary(mod.plot)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm * log.lifespan 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   235.0328 251.9059 -111.5164
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.5521672 
> 
> Coefficients:
>                                     Value Std.Error   t-value p-value
> (Intercept)                      3.222911 0.4660986  6.914656  0.0000
> Endo_ectothermendo              -0.103630 0.5207431 -0.199003  0.8426
> log.lifespan                    -0.086519 0.1295737 -0.667717  0.5056
> Endo_ectothermendo:log.lifespan -0.234612 0.1649920 -1.421960  0.1577
> 
>  Correlation: 
>                                 (Intr) End_ct lg.lfs
> Endo_ectothermendo              -0.580              
> log.lifespan                    -0.798  0.676       
> Endo_ectothermendo:log.lifespan  0.619 -0.912 -0.778
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -2.41100472  0.05857414  0.58071564  1.08908175  2.41709030 
> 
> Residual standard error: 0.7673209 
> Degrees of freedom: 123 total; 119 residual
## Model diagnosis


layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))



## Checking homogeneity of variance


plot(fitted(pagel.mod.red), resid(pagel.mod.red), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()

par(new = TRUE)
plot(fitted(pagel.mod.red), resid(pagel.mod.red), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)

abline(h = 0, col = "darkorange", lwd = 2)


## Checking normality

qqnorm(resid(pagel.mod.red), col = "darkgrey", type = 'n', las = 1)
grid()

par(new = TRUE)
qqnorm(resid(pagel.mod.red), col = "darkgrey", las = 1)
qqline(resid(pagel.mod.red), col = "dodgerblue", lwd = 2)



Figure 3

#png("figure3-revised.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure3-revised.pdf")

plot(log(Average_Telomere_Length_kb[Endo_ectotherm == 'ecto']) ~ log.lifespan[Endo_ectotherm == 'ecto'], data = pruned_dat, pch = 21, bg = alpha("purple", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)", type = 'n')

grid(nx = NULL, ny = NULL, col = "lightgray", lwd = 1)
par(new = TRUE)

plot(log(Average_Telomere_Length_kb[Endo_ectotherm == 'ecto']) ~ log.lifespan[Endo_ectotherm == 'ecto'], data = pruned_dat, pch = 21, bg = alpha("purple", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)")
points(log(Average_Telomere_Length_kb[Endo_ectotherm == 'endo']) ~ log.lifespan[Endo_ectotherm == 'endo'], data = pruned_dat, pch = 21, bg = alpha("orange", 0.5), las = 1, xlab = "Lifespan (log yrs)", ylab = "Telomere length (log kb)")

SSX <- sum(round((log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']) - mean(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto'])))^2), 2)
s2 <- var(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']))
n <- length(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']))
x <- seq(min(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'ecto']), max(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'ecto']), 0.06)
m.x <- mean(round(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'ecto']), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(mod.plot)[1] + coef(mod.plot)[2]*x) + ic.s
lower.i <- (coef(mod.plot)[1] + coef(mod.plot)[2]*x) + ic.i

##par(new = TRUE)

polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("purple", 0.3))
lines(x = x, y = (coef(mod.plot)[1] + (coef(mod.plot)[2]*x)), lwd = 2, col = alpha("purple", 0.5))

SSX <- sum(round((log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']) - mean(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo'])))^2), 2)
s2 <- var(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']))
n <- length(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']))
x <- seq(min(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'endo']), max(pruned_dat$log.lifespan[pruned_dat$Endo_ectotherm == 'endo']), 0.06)
m.x <- mean(round(log(pruned_dat$Average_Telomere_Length_kb[pruned_dat$Endo_ectotherm == 'endo']), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x) + ic.s
lower.i <- ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x) + ic.i

polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("orange", 0.3))
lines(x = x, y = ((coef(mod.plot)[1] + coef(mod.plot)[3]) + (coef(mod.plot)[2] + coef(mod.plot)[4])*x), lwd = 2, col = alpha("orange", 0.5))
##summary(mod.plot)

legend('bottomleft', legend = c('Ectotherms', 'Endotherms'), pch = 16, col = c('purple', 'orange'), lwd = 1, bty = 'n', cex = 0.8)

#dev.off()

Figure 2

## Phylogenetic signal (Blomberg et al.'s K statitics

telo.length <- setNames(log(pruned_dat$Average_Telomere_Length_kb), rownames(pruned_dat))
telo.length
>       Dicentrarchus_labrax       Pseudocaranx_wrighti       Xiphophorus_hellerii      Xiphophorus_maculatus          Carassius_auratus 
>                  1.2354715                  1.4011830                  1.4469190                  1.4469190                  1.5260563 
>            Cyprinus_carpio  Lutjanus_argentimaculatus     Upeneichthys_vlamingii   Acanthopagrus_schlegelii          Macquaria_ambigua 
>                  1.5303947                  1.6601310                  1.7011051                  1.7155981                  1.8148247 
>     Nothobranchius_furzeri    Nothobranchius_rachovii               Gadus_morhua    Platycephalus_bassensis                Danio_rerio 
>                  1.8885837                  2.0502702                  2.4313300                  2.5006159                  2.7725887 
>         Alligator_sinensis             Lacerta_agilis           Zootoca_vivipara           Emys_orbicularis        Oncorhynchus_mykiss 
>                  2.7825391                  2.9231616                  2.9606231                  2.9957323                  2.9957323 
>      Neoceratodus_forsteri   Urolophus_paucimaculatus        Callorhinchus_milii              Liasis_fuscus      Meriones_unguiculatus 
>                  3.1871787                  3.1941732                  3.2023399                  3.2976870                  3.3672958 
> Alligator_mississippiensis        Chinchilla_lanigera         Ondatra_zibethicus       Mesocricetus_auratus          Trachemys_scripta 
>                  3.4307562                  3.6375862                  3.8066625                  3.9120230                  3.9120230 
>           Myocastor_coypus              Marmota_monax            Tamias_striatus         Accipiter_gentilis   Haliaeetus_leucocephalus 
>                  3.9889840                  4.0943446                  4.1431347                  0.2623643                  0.7419373 
>        Chrysolophus_pictus     Zalophus_californianus    Aphelocoma_coerulescens          Castor_canadensis           Ateles_geoffroyi 
>                  1.1631508                  1.6094379                  1.7101878                  1.9459101                  1.9459101 
>              Fregata_minor         Pygoscelis_adeliae                Uria_lomvia            Crocuta_crocuta      Macronectes_giganteus 
>                  1.9586853                  1.9947003                  2.0412203                  2.0794415                  2.0893919 
>          Macronectes_halli       Calonectris_diomedea         Balaena_mysticetus           Saimiri_sciureus       Pteropus_rodricensis 
>                  2.1488510                  2.1587147                  2.1972246                  2.1972246                  2.1972246 
>               Homo_sapiens            Aplodontia_rufa     Peromyscus_maniculatus      Tachycineta_albilinea            Hirundo_rustica 
>                  2.1972246                  2.1972246                  2.1972246                  2.2016592                  2.2586332 
>      Haematopus_ostralegus           Rissa_tridactyla           Diomedea_exulans               Pan_paniscus  Hydrochoerus_hydrochaeris 
>                  2.2813615                  2.2828925                  2.2925348                  2.3025851                  2.3025851 
>     Giraffa_camelopardalis             Pongo_pygmaeus        Ceratotherium_simum            Vireo_olivaceus              Turdus_merula 
>                  2.3025851                  2.3025851                  2.3025851                  2.3841651                  2.4024304 
>        Larus_crassirostris            Pan_troglodytes  Passerculus_sandwichensis           Branta_leucopsis          Amazona_amazonica 
>                  2.4362415                  2.4510051                  2.4535880                  2.4638532                  2.4680995 
>    Myrmecophaga_tridactyla            Tapirus_indicus            Ursus_maritimus               Equus_grevyi             Sterna_hirundo 
>                  2.4849066                  2.4849066                  2.4849066                  2.4849066                  2.5494452 
>      Eschrichtius_robustus  Chaetophractus_vellerosus          Buteo_jamaicensis     Aphelocoma_ultramarina        Tachycineta_bicolor 
>                  2.5649494                  2.5649494                  2.5877640                  2.6100698                  2.6100698 
>        Taeniopygia_guttata            Calidris_alpina        Meleagris_gallopavo         Loxodonta_africana          Muntiacus_reevesi 
>                  2.6173958                  2.6253930                  2.6318888                  2.6390573                  2.6390573 
>          Muntiacus_muntjak         Fulmarus_glacialis             Cepphus_grylle                  Sula_sula          Macaca_nemestrina 
>                  2.6390573                  2.6567569                  2.6803364                  2.6946272                  2.7013612 
>            Elephas_maximus                Tupaia_tana          Procavia_capensis        Macaca_fascicularis            Gorilla_gorilla 
>                  2.7080502                  2.7080502                  2.7080502                  2.7146947                  2.7343675 
>        Strigops_habroptila      Heterocephalus_glaber             Macaca_mulatta           Lonchura_striata         Tursiops_truncatus 
>                  2.7507904                  2.7725887                  2.7725887                  2.7948393                  2.8332133 
>        Choloepus_hoffmanni         Tachymarptis_melba                Lemur_catta      Oceanodroma_leucorhoa          Ochotona_princeps 
>                  2.8332133                  2.9329522                  2.9444390                  2.9460167                  2.9957323 
>            Riparia_riparia        Colinus_virginianus            Ailurus_fulgens      Tadarida_brasiliensis           Myotis_lucifugus 
>                  3.0568274                  3.1527360                  3.2188758                  3.2580965                  3.4011974 
>           Falco_sparverius       Didelphis_virginiana       Galegeeska_rufescens Macroscelides_proboscideus       Atelerix_albiventris 
>                  3.4078419                  3.5553481                  3.6109179                  3.6109179                  3.6375862 
>            Buteo_swainsoni                Parus_major         Lepus_californicus       Sciurus_carolinensis      Oryctolagus_cuniculus 
>                  3.7471484                  3.8971124                  3.9120230                  3.9120230                  3.9120230 
>            Setifer_setosus       Sylvilagus_aquaticus            Panthera_tigris 
>                  3.9120230                  3.9120230                  3.9120230
k_tl <- phylosig(pruned_data_tree, telo.length, test = TRUE, nsim = 10000)
attributes(k_tl)
> $names
> [1] "K"     "P"     "sim.K"
> 
> $class
> [1] "phylosig"
> 
> $method
> [1] "K"
> 
> $test
> [1] TRUE
> 
> $se
> [1] FALSE
head(k_tl$sim.K)
> [1] 0.16697908 0.08027691 0.10675328 0.05298757 0.05080825 0.03556077
k_tl$K
> [1] 0.1669791
k_tl$P
> [1] 2e-04
#png("figure2.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure2.pdf")


plot(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "", type = 'n', axes = FALSE)
grid()

par(new = TRUE)
hist(k_tl$sim.K, bty = "o", ylim = c(0, 3000), las = 1, ylab = "Null distribution of K", xlab = "K", main = "")

abline(v = k_tl$K, lwd = 2, lty = "dotted")
text(x = k_tl$K - 0.02, y = 2500, "Observed value \n of K")
box()

#dev.off()

## Estimating lambda from PGLS models

pagel.model <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0.8, phy = pruned_data_tree, form = ~Scientific_name), data = pruned_dat, method = "ML")
summary(pagel.model)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   235.1572 252.0303 -111.5786
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
>    lambda 
> 0.6254488 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.629715 0.3916907  9.266788  0.0000
> Endo_ectothermendo -0.819368 0.2207192 -3.712263  0.0003
> log.lifespan       -0.134538 0.1130636 -1.189932  0.2364
> log.mass           -0.033764 0.0258447 -1.306413  0.1939
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.041              
> log.lifespan       -0.502 -0.101       
> log.mass            0.075  0.013 -0.678
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -2.34917753 -0.02875553  0.50192730  0.92152007  2.39155284 
> 
> Residual standard error: 0.8117616 
> Degrees of freedom: 123 total; 119 residual
intervals(pagel.model, which = 'var-cov')$corStruct
>            lower      est.     upper
> lambda 0.3607961 0.6254488 0.8901015
> attr(,"label")
> [1] "Correlation structure:"
pagel.model0 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   281.5872 295.6481 -135.7936
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>      0 
> 
> Coefficients:
>                         Value  Std.Error   t-value p-value
> (Intercept)         3.0127181 0.27461330 10.970765  0.0000
> Endo_ectothermendo -0.0099686 0.15910482 -0.062654  0.9501
> log.lifespan       -0.1146562 0.11583847 -0.989794  0.3243
> log.mass           -0.0007418 0.02556906 -0.029011  0.9769
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.108              
> log.lifespan       -0.748 -0.291       
> log.mass            0.131  0.098 -0.650
> 
> Standardized residuals:
>        Min         Q1        Med         Q3        Max 
> -3.2551171 -0.5625105 -0.0212407  0.5130302  1.9228926 
> 
> Residual standard error: 0.7298431 
> Degrees of freedom: 123 total; 119 residual
intervals(pagel.model0, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>     lower      est.     upper 
> 0.6489217 0.7298431 0.8340042
pagel.model0.5 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(0.5, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model0.5)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   233.9806 248.0415 -111.9903
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>    0.5 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.566355 0.3532026 10.097194  0.0000
> Endo_ectothermendo -0.752855 0.2101979 -3.581649  0.0005
> log.lifespan       -0.131519 0.1105688 -1.189477  0.2366
> log.mass           -0.031862 0.0254440 -1.252233  0.2129
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.053              
> log.lifespan       -0.542 -0.101       
> log.mass            0.094  0.003 -0.685
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -2.57822131 -0.06855442  0.51155121  0.98428268  2.56306826 
> 
> Residual standard error: 0.7436376 
> Degrees of freedom: 123 total; 119 residual
intervals(pagel.model0.5, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>     lower      est.     upper 
> 0.6611868 0.7436376 0.8497675
pagel.model1 <- gls(log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan + log.mass, correlation = corPagel(1, phy = pruned_data_tree, form = ~Scientific_name, fixed = TRUE), data = pruned_dat, method = "ML")
summary(pagel.model1)
> Generalized least squares fit by maximum likelihood
>   Model: log(Average_Telomere_Length_kb) ~ Endo_ectotherm + log.lifespan +      log.mass 
>   Data: pruned_dat 
>        AIC      BIC    logLik
>   272.3855 286.4464 -131.1927
> 
> Correlation Structure: corPagel
>  Formula: ~Scientific_name 
>  Parameter estimate(s):
> lambda 
>      1 
> 
> Coefficients:
>                        Value Std.Error   t-value p-value
> (Intercept)         3.901248 0.8933978  4.366753  0.0000
> Endo_ectothermendo -1.172484 0.2869553 -4.085949  0.0001
> log.lifespan       -0.142480 0.1268266 -1.123421  0.2635
> log.mass           -0.048054 0.0328901 -1.461046  0.1466
> 
>  Correlation: 
>                    (Intr) End_ct lg.lfs
> Endo_ectothermendo -0.025              
> log.lifespan       -0.268 -0.136       
> log.mass           -0.086  0.108 -0.499
> 
> Standardized residuals:
>         Min          Q1         Med          Q3         Max 
> -1.11576421  0.04130495  0.29665511  0.48025508  1.15040877 
> 
> Residual standard error: 1.926538 
> Degrees of freedom: 123 total; 119 residual
intervals(pagel.model1, which = 'var-cov')
> Approximate 95% confidence intervals
> 
>  Residual standard error:
>    lower     est.    upper 
> 1.712933 1.926538 2.201488
anova(pagel.model, pagel.model0)
>              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> pagel.model      1  6 235.1572 252.0303 -111.5786                        
> pagel.model0     2  5 281.5872 295.6481 -135.7936 1 vs 2 48.42996  <.0001
anova(pagel.model, pagel.model0.5)
>                Model df      AIC      BIC    logLik   Test   L.Ratio p-value
> pagel.model        1  6 235.1572 252.0303 -111.5786                         
> pagel.model0.5     2  5 233.9806 248.0415 -111.9903 1 vs 2 0.8233961  0.3642
anova(pagel.model0, pagel.model0.5)
>                Model df      AIC      BIC    logLik
> pagel.model0       1  5 281.5872 295.6481 -135.7936
> pagel.model0.5     2  5 233.9806 248.0415 -111.9903
anova(pagel.model, pagel.model1)
>              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> pagel.model      1  6 235.1572 252.0303 -111.5786                        
> pagel.model1     2  5 272.3855 286.4464 -131.1927 1 vs 2 39.22828  <.0001
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 0] <- "absent"
pruned_dat$Telomerase_activity[pruned_dat$Telomerase_activity == 1] <- "present"
pruned_dat$Telomerase_activity[is.na(pruned_dat$Telomerase_activity)] <- "N/A"

tel.act <- setNames(pruned_dat$Telomerase_activity, rownames(pruned_dat))
tel.act
>       Dicentrarchus_labrax       Pseudocaranx_wrighti       Xiphophorus_hellerii      Xiphophorus_maculatus          Carassius_auratus 
>                      "N/A"                      "N/A"                      "N/A"                      "N/A"                      "N/A" 
>            Cyprinus_carpio  Lutjanus_argentimaculatus     Upeneichthys_vlamingii   Acanthopagrus_schlegelii          Macquaria_ambigua 
>                      "N/A"                      "N/A"                      "N/A"                      "N/A"                      "N/A" 
>     Nothobranchius_furzeri    Nothobranchius_rachovii               Gadus_morhua    Platycephalus_bassensis                Danio_rerio 
>                  "present"                      "N/A"                  "present"                      "N/A"                  "present" 
>         Alligator_sinensis             Lacerta_agilis           Zootoca_vivipara           Emys_orbicularis        Oncorhynchus_mykiss 
>                   "absent"                      "N/A"                      "N/A"                      "N/A"                  "present" 
>      Neoceratodus_forsteri   Urolophus_paucimaculatus        Callorhinchus_milii              Liasis_fuscus      Meriones_unguiculatus 
>                      "N/A"                      "N/A"                      "N/A"                      "N/A"                  "present" 
> Alligator_mississippiensis        Chinchilla_lanigera         Ondatra_zibethicus       Mesocricetus_auratus          Trachemys_scripta 
>                      "N/A"                  "present"                  "present"                  "present"                      "N/A" 
>           Myocastor_coypus              Marmota_monax            Tamias_striatus         Accipiter_gentilis   Haliaeetus_leucocephalus 
>                  "present"                  "present"                  "present"                      "N/A"                      "N/A" 
>        Chrysolophus_pictus     Zalophus_californianus    Aphelocoma_coerulescens          Castor_canadensis           Ateles_geoffroyi 
>                      "N/A"                   "absent"                      "N/A"                   "absent"                   "absent" 
>              Fregata_minor         Pygoscelis_adeliae                Uria_lomvia            Crocuta_crocuta      Macronectes_giganteus 
>                      "N/A"                      "N/A"                      "N/A"                   "absent"                      "N/A" 
>          Macronectes_halli       Calonectris_diomedea         Balaena_mysticetus           Saimiri_sciureus       Pteropus_rodricensis 
>                      "N/A"                      "N/A"                   "absent"                   "absent"                   "absent" 
>               Homo_sapiens            Aplodontia_rufa     Peromyscus_maniculatus      Tachycineta_albilinea            Hirundo_rustica 
>                   "absent"                   "absent"                  "present"                      "N/A"                      "N/A" 
>      Haematopus_ostralegus           Rissa_tridactyla           Diomedea_exulans               Pan_paniscus  Hydrochoerus_hydrochaeris 
>                      "N/A"                      "N/A"                      "N/A"                   "absent"                   "absent" 
>     Giraffa_camelopardalis             Pongo_pygmaeus        Ceratotherium_simum            Vireo_olivaceus              Turdus_merula 
>                   "absent"                   "absent"                   "absent"                      "N/A"                      "N/A" 
>        Larus_crassirostris            Pan_troglodytes  Passerculus_sandwichensis           Branta_leucopsis          Amazona_amazonica 
>                      "N/A"                      "N/A"                      "N/A"                      "N/A"                      "N/A" 
>    Myrmecophaga_tridactyla            Tapirus_indicus            Ursus_maritimus               Equus_grevyi             Sterna_hirundo 
>                   "absent"                   "absent"                   "absent"                   "absent"                  "present" 
>      Eschrichtius_robustus  Chaetophractus_vellerosus          Buteo_jamaicensis     Aphelocoma_ultramarina        Tachycineta_bicolor 
>                   "absent"                   "absent"                      "N/A"                      "N/A"                  "present" 
>        Taeniopygia_guttata            Calidris_alpina        Meleagris_gallopavo         Loxodonta_africana          Muntiacus_reevesi 
>                  "present"                      "N/A"                      "N/A"                   "absent"                   "absent" 
>          Muntiacus_muntjak         Fulmarus_glacialis             Cepphus_grylle                  Sula_sula          Macaca_nemestrina 
>                   "absent"                      "N/A"                      "N/A"                      "N/A"                      "N/A" 
>            Elephas_maximus                Tupaia_tana          Procavia_capensis        Macaca_fascicularis            Gorilla_gorilla 
>                   "absent"                  "present"                   "absent"                   "absent"                      "N/A" 
>        Strigops_habroptila      Heterocephalus_glaber             Macaca_mulatta           Lonchura_striata         Tursiops_truncatus 
>                      "N/A"                  "present"                   "absent"                      "N/A"                   "absent" 
>        Choloepus_hoffmanni         Tachymarptis_melba                Lemur_catta      Oceanodroma_leucorhoa          Ochotona_princeps 
>                   "absent"                      "N/A"                   "absent"                  "present"                  "present" 
>            Riparia_riparia        Colinus_virginianus            Ailurus_fulgens      Tadarida_brasiliensis           Myotis_lucifugus 
>                      "N/A"                      "N/A"                   "absent"                  "present"                  "present" 
>           Falco_sparverius       Didelphis_virginiana       Galegeeska_rufescens Macroscelides_proboscideus       Atelerix_albiventris 
>                      "N/A"                  "present"                  "present"                  "present"                  "present" 
>            Buteo_swainsoni                Parus_major         Lepus_californicus       Sciurus_carolinensis      Oryctolagus_cuniculus 
>                      "N/A"                      "N/A"                   "absent"                  "present"                   "absent" 
>            Setifer_setosus       Sylvilagus_aquaticus            Panthera_tigris 
>                  "present"                   "absent"                  "present"
TA <- to.matrix(tel.act, unique(tel.act))
TA <- TA[pruned_data_tree$tip.label, ]

life.span <- setNames(pruned_dat$log.lifespan, rownames(pruned_dat))
log_mass <- setNames(pruned_dat$log.mass, rownames(pruned_dat))
telo.length <- setNames(log(pruned_dat$Average_Telomere_Length_kb), rownames(pruned_dat))
plotTree(pruned_data_tree, ftype = "off", mar = c(3, 2, 2, 3))

tiplabels(pie = TA, piecol = c("white", "gray", "black"), cex = 0.22, offset = 4.3)

par(xpd = TRUE)

legend(x = 150, y = 0.2, legend = unique(tel.act), pch = 21, pt.bg = c("white", "gray", "black"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)


par(new = TRUE)
par(mar = c(3, 32, 2, 1.1))
barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)

## Ancestral state reconstruction of telomere size

fit <- fastAnc(pruned_data_tree, telo.length, vars = TRUE, CI = TRUE)
fit$CI[1, ]
> [1] 1.156184 4.676634
obj <- contMap(pruned_data_tree, telo.length, plot = FALSE)



Figure 1

##png("figure1.png", width = 7, height = 7, units = "in", res = 360)
##pdf("figure1.pdf")


plot(obj, ftype = "off", legend = FALSE, ylim = c(1-0.09*(Ntip(obj$tree)-1), Ntip(obj$tree)), mar = c(1, 0.1, 1, 4), lwd = 1.5)
add.color.bar(150, obj$cols,title = "Log telomere length (kb)", lims = obj$lims, digits = 3, prompt = FALSE, x = 0,
              y = 1-0.08*(Ntip(obj$tree)-1), lwd = 4, fsize = 0.6, subtitle = "")

par(xpd = TRUE)

text(300, 119, 'Testudines')
text(300, 101, 'Aves', col = 'red')
text(300, 82, 'Crocodilia')
text(240, 74, 'Squamata')
text(240, 40, 'Mammalia', col = 'red')
text(160, 15, 'Osteichthyes')
text(160, 4.5, 'Chondrichthyes')


#cladelabels(pruned_data_tree, node = 217, "Ectotherms", offset = 6)
#cladelabels(pruned_data_tree, node = 153, "Endotherms", offset = 6)
#segments(518, 0, 518, 18)
#segments(518, 0, 508, 0)
#segments(518, 18, 508, 18)
#text(528, 9, 'Ectotherms', srt = 90)


tiplabels(pie = TA, piecol = c("gray", "black", "white"), cex = 0.16, offset = 5.7)
legend(x = 200, y = -5, legend = unique(tel.act), pch = 21, pt.bg = c("gray", "black", "white"), pt.cex = 1, bty = "n", title = "Telomerase activity", cex = 0.7, horiz = TRUE)

par(new = TRUE)
par(mar = c(3.6, 30.8, 3, 2.2))

barplot(life.span[pruned_data_tree$tip.label], horiz = TRUE, width = 1.07, space = 0,
        ylim = c(1, length(pruned_data_tree$tip.label))-0.5, names = "", las = 2, cex.axis = 0.5, axes = FALSE)

axis(1, at = round(seq(min(life.span), max(life.span), 1.5), 1), labels = FALSE)
text(round(seq(min(life.span), max(life.span), 1.5), 1), par("usr")[3] - 0.2, labels = round(seq(min(life.span), max(life.span), 1.5), 1), srt = 50, pos = 1, xpd = TRUE, cex = 0.5, offset = 1)
mtext("Log lifespan \n (years)", side = 1, line = 1.6, cex = 0.5, font = 2)

##dev.off()
## Correlated evolution under the threshold model


## Removing NAs from the dataset

pruned_dat$Telomerase_activity
>   [1] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "present" "N/A"     "present"
>  [14] "N/A"     "present" "absent"  "N/A"     "N/A"     "N/A"     "present" "N/A"     "N/A"     "N/A"     "N/A"     "present" "N/A"    
>  [27] "present" "present" "present" "N/A"     "present" "present" "present" "N/A"     "N/A"     "N/A"     "absent"  "N/A"     "absent" 
>  [40] "absent"  "N/A"     "N/A"     "N/A"     "absent"  "N/A"     "N/A"     "N/A"     "absent"  "absent"  "absent"  "absent"  "absent" 
>  [53] "present" "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "absent"  "absent"  "absent"  "absent"  "absent"  "N/A"     "N/A"    
>  [66] "N/A"     "N/A"     "N/A"     "N/A"     "N/A"     "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "N/A"    
>  [79] "N/A"     "present" "present" "N/A"     "N/A"     "absent"  "absent"  "absent"  "N/A"     "N/A"     "N/A"     "N/A"     "absent" 
>  [92] "present" "absent"  "absent"  "N/A"     "N/A"     "present" "absent"  "N/A"     "absent"  "absent"  "N/A"     "absent"  "present"
> [105] "present" "N/A"     "N/A"     "absent"  "present" "present" "N/A"     "present" "present" "present" "present" "N/A"     "N/A"    
> [118] "absent"  "present" "absent"  "present" "absent"  "present"
pruned_dat_not_NAs <- pruned_dat[!pruned_dat$Telomerase_activity == "N/A", ]
pruned_dat_not_NAs$Telomerase_activity
>  [1] "present" "present" "present" "absent"  "present" "present" "present" "present" "present" "present" "present" "present" "absent" 
> [14] "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "absent"  "absent" 
> [27] "absent"  "absent"  "absent"  "absent"  "absent"  "present" "absent"  "absent"  "present" "present" "absent"  "absent"  "absent" 
> [40] "absent"  "present" "absent"  "absent"  "present" "absent"  "absent"  "absent"  "absent"  "present" "present" "absent"  "present"
> [53] "present" "present" "present" "present" "present" "absent"  "present" "absent"  "present" "absent"  "present"
check2 <- name.check(pruned_data_tree, pruned_dat_not_NAs)
rm_phy2 <- check2$tree_not_data
rm_dat2 <- check2$data_not_tree
pruned_data_tree_not_NAs <- drop.tip(pruned_data_tree, rm_phy2)
pruned_dat_not_NAs <- subset(pruned_dat_not_NAs, subset = pruned_dat_not_NAs$Scientific_name %in% pruned_data_tree_not_NAs$tip, select = names(pruned_dat_not_NAs))

pruned_dat_not_NAs$Telomerase_activity <- as.factor(pruned_dat_not_NAs$Telomerase_activity)
str(pruned_dat_not_NAs)
> 'data.frame': 63 obs. of  29 variables:
>  $ Common.name               : chr  "African killifish" "Atlantic cod" "Zebrafish" "Chinese alligator" ...
>  $ Scientific_name           : chr  "Nothobranchius_furzeri" "Gadus_morhua" "Danio_rerio" "Alligator_sinensis" ...
>  $ Domain                    : chr  "Eukaryota" "Eukaryota" "Eukaryota" "Eukaryota" ...
>  $ Kingdom                   : chr  "Metazoa" "Metazoa" "Metazoa" "Metazoa" ...
>  $ Phylum                    : chr  "Chordata" "Chordata" "Chordata" "Chordata" ...
>  $ Class                     : chr  "Fish" "Fish" "Fish" "Reptilia" ...
>  $ Order                     : chr  "Cyprinodontiformes" "Gadiformes" "Cypriniformes" "Crocodilia" ...
>  $ Family                    : chr  "Nothobranchiidae" "Gadidae" "Cyprinidae" "Alligatoridae" ...
>  $ Genus                     : chr  "Nothobranchius" "Gadus" "Danio" "Alligator" ...
>  $ Species                   : chr  "furzeri" "morhua" "rerio" "sinensis" ...
>  $ Endo_ectotherm            : chr  "ecto" "ecto" "ecto" "ecto" ...
>  $ Adult_mass_grams          : num  2 52800 0.5 14600 12235 ...
>  $ log_mass                  : num  0.301 4.723 -0.301 4.164 4.088 ...
>  $ Lifespan_years            : num  0.23 25 5.5 65 11 6.3 17.2 10 3.9 19 ...
>  $ Loglifespan               : num  -0.638 1.398 0.74 1.813 1.041 ...
>  $ Average_Telomere_Length_kb: num  6.61 11.37 16 16.16 20 ...
>  $ Telomerase_activity       : Factor w/ 2 levels "absent","present": 2 2 2 1 2 2 2 2 2 2 ...
>  $ Tissue.type.for.TA        : chr  "Various tissues" "Multiple somatic tissues" "Various tissues" "Blood" ...
>  $ Source.for.TA             : chr  "Hartmann et al. 2009" "Abechuco et al. 2014" "Achelin et al. 2011" "Guo et al. 2023" ...
>  $ Tissue.type.for.TL        : chr  "Muscle and skin" "Multiple somatic tissues" "Various tissues" "Blood" ...
>  $ Methodology.for.TL        : chr  "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" "TRF (denatured)" ...
>  $ TRF_type                  : chr  "A" "A" "A" "A" ...
>  $ Source.for.TL             : chr  "Hartmann et al. 2009" "de Abechuco et al. 2016" "Lund et al. 2009; Achelin et al. 2011" "Xu et al. 2009" ...
>  $ Source.for.Lifespan       : chr  "Valdesalici and Cellerino 2003" "AnAge" "AnAge" "AnAge" ...
>  $ Source.for.Mass           : chr  "Pola?ik and Reichard 2010" "AnAge" "Fishbase" "AnAge" ...
>  $ Captivity                 : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Notes                     : chr  "" "" "" "" ...
>  $ log.lifespan              : num  0.207 3.258 1.872 4.19 2.485 ...
>  $ log.mass                  : num  1.099 10.874 0.405 9.589 9.412 ...
head(pruned_dat_not_NAs)
>                              Common.name        Scientific_name    Domain Kingdom   Phylum    Class              Order           Family
> Nothobranchius_furzeri African killifish Nothobranchius_furzeri Eukaryota Metazoa Chordata     Fish Cyprinodontiformes Nothobranchiidae
> Gadus_morhua                Atlantic cod           Gadus_morhua Eukaryota Metazoa Chordata     Fish         Gadiformes          Gadidae
> Danio_rerio                    Zebrafish            Danio_rerio Eukaryota Metazoa Chordata     Fish      Cypriniformes       Cyprinidae
> Alligator_sinensis     Chinese alligator     Alligator_sinensis Eukaryota Metazoa Chordata Reptilia         Crocodilia    Alligatoridae
> Oncorhynchus_mykiss        Rainbow trout    Oncorhynchus_mykiss Eukaryota Metazoa Chordata     Fish      Salmoniformes       Salmonidae
> Meriones_unguiculatus   Mongolian Gerbil  Meriones_unguiculatus Eukaryota Metazoa Chordata Mammalia           Rodentia          Muridae
>                                 Genus      Species Endo_ectotherm Adult_mass_grams  log_mass Lifespan_years Loglifespan
> Nothobranchius_furzeri Nothobranchius      furzeri           ecto              2.0  0.301030           0.23  -0.6382722
> Gadus_morhua                    Gadus       morhua           ecto          52800.0  4.722634          25.00   1.3979400
> Danio_rerio                     Danio        rerio           ecto              0.5 -0.301030           5.50   0.7403627
> Alligator_sinensis          Alligator     sinensis           ecto          14600.0  4.164353          65.00   1.8129134
> Oncorhynchus_mykiss      Oncorhynchus       mykiss           ecto          12235.0  4.087604          11.00   1.0413927
> Meriones_unguiculatus        Meriones unguiculatus           ecto             53.2  1.725912           6.30   0.7993405
>                        Average_Telomere_Length_kb Telomerase_activity       Tissue.type.for.TA        Source.for.TA
> Nothobranchius_furzeri                      6.610             present          Various tissues Hartmann et al. 2009
> Gadus_morhua                               11.374             present Multiple somatic tissues Abechuco et al. 2014
> Danio_rerio                                16.000             present          Various tissues  Achelin et al. 2011
> Alligator_sinensis                         16.160              absent                    Blood      Guo et al. 2023
> Oncorhynchus_mykiss                        20.000             present         Multiple tissues Panasiak et al. 2023
> Meriones_unguiculatus                      29.000             present Multiple somatic tissues Seluanov et al. 2007
>                              Tissue.type.for.TL Methodology.for.TL TRF_type                         Source.for.TL
> Nothobranchius_furzeri          Muscle and skin    TRF (denatured)        A                  Hartmann et al. 2009
> Gadus_morhua           Multiple somatic tissues    TRF (denatured)        A               de Abechuco et al. 2016
> Danio_rerio                     Various tissues    TRF (denatured)        A Lund et al. 2009; Achelin et al. 2011
> Alligator_sinensis                        Blood    TRF (denatured)        A                        Xu et al. 2009
> Oncorhynchus_mykiss                       Blood    TRF (denatured)        A                   Lejnine et al. 1995
> Meriones_unguiculatus  Multiple somatic tissues    TRF (denatured)        A                  Seluanov et al. 2007
>                                   Source.for.Lifespan           Source.for.Mass Captivity Notes log.lifespan   log.mass
> Nothobranchius_furzeri Valdesalici and Cellerino 2003 Pola?ik and Reichard 2010         0          0.2070142  1.0986123
> Gadus_morhua                                    AnAge                     AnAge         0          3.2580965 10.8742854
> Danio_rerio                                     AnAge                  Fishbase         0          1.8718022  0.4054651
> Alligator_sinensis                              AnAge                     AnAge         0          4.1896547  9.5888453
> Oncorhynchus_mykiss                             AnAge                     AnAge         0          2.4849066  9.4121377
> Meriones_unguiculatus                           AnAge                     AnAge         0          1.9878743  3.9926809
pruned_data_tree_not_NAs
> 
> Phylogenetic tree with 63 tips and 62 internal nodes.
> 
> Tip labels:
>   Oncorhynchus_mykiss, Nothobranchius_furzeri, Gadus_morhua, Danio_rerio, Didelphis_virginiana, Atelerix_albiventris, ...
> Node labels:
>   '12', '48', '57', '49', '194', '166', ...
> 
> Rooted; includes branch lengths.
name.check(pruned_data_tree_not_NAs, pruned_dat_not_NAs)
> [1] "OK"
names(pruned_dat_not_NAs)
>  [1] "Common.name"                "Scientific_name"            "Domain"                     "Kingdom"                   
>  [5] "Phylum"                     "Class"                      "Order"                      "Family"                    
>  [9] "Genus"                      "Species"                    "Endo_ectotherm"             "Adult_mass_grams"          
> [13] "log_mass"                   "Lifespan_years"             "Loglifespan"                "Average_Telomere_Length_kb"
> [17] "Telomerase_activity"        "Tissue.type.for.TA"         "Source.for.TA"              "Tissue.type.for.TL"        
> [21] "Methodology.for.TL"         "TRF_type"                   "Source.for.TL"              "Source.for.Lifespan"       
> [25] "Source.for.Mass"            "Captivity"                  "Notes"                      "log.lifespan"              
> [29] "log.mass"
## Set the number of generations

ngen <- 5e6

## Run MCMC

pruned_dat_not_NAs$Average_Telomere_Length_kb <- log(pruned_dat_not_NAs$Average_Telomere_Length_kb)

mcmc.model <- threshBayes(pruned_data_tree_not_NAs, pruned_dat_not_NAs[ , c(16, 17)], type = c("cont", "disc"), ngen = ngen,
                          plot = FALSE, control = list(print.interval = 5e+05))
> Starting MCMC....
> generation: 500000; mean acceptance rate: 0.23
> generation: 1000000; mean acceptance rate: 0.24
> generation: 1500000; mean acceptance rate: 0.27
> generation: 2000000; mean acceptance rate: 0.2
> generation: 2500000; mean acceptance rate: 0.26
> generation: 3000000; mean acceptance rate: 0.23
> generation: 3500000; mean acceptance rate: 0.21
> generation: 4000000; mean acceptance rate: 0.23
> generation: 4500000; mean acceptance rate: 0.21
> generation: 5000000; mean acceptance rate: 0.22
> Done MCMC.
mcmc.model
> 
> Object of class "threshBayes" consisting of a matrix (L) of
> sampled liabilities for the tips of the tree & a second matrix
> (par) with the sample model parameters & correlation.
> 
> Mean correlation (r) from the posterior sample is: 0.61357.
> 
> Ordination of discrete traits:
> 
>   Trait 2: absent <-> present
## Pull out the post burn-in sample and compute HPD

r.mcmc <- tail(mcmc.model$par$r, 0.8*nrow(mcmc.model$par))
class(r.mcmc) <- "mcmc"

hpd.r <- HPDinterval(r.mcmc)
hpd.r
>          lower     upper
> var1 0.3487216 0.8541938
> attr(,"Probability")
> [1] 0.9500012



Figure S3

## Profile plots from a Bayesian MCMC analysis of the threshold model

#png("figureS3.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figureS3.pdf")



layout(matrix(c(0, 0, 0, 0,
                1, 1, 1, 1,
                1, 1, 1, 1,
                2, 2, 2, 2,
                2, 2, 2, 2,
                0, 0, 0, 0), nrow = 6, ncol = 4, byrow = TRUE))

plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'n', ylab = 'Log(L)', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$logL ~ mcmc.model$par$gen, col = 'gray', type = 'l', ylab = 'Log(L)', xlab = 'Generation', las = 1)
mtext('(a)', side = 2, line = 1.5, at = -160, las = 1)

plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'n', ylab = 'r', xlab = 'Generation', las = 1)
grid()
par(new = TRUE)
plot(mcmc.model$par$r ~ mcmc.model$par$gen, col = 'purple', type = 'l', ylab = 'r', xlab = 'Generation', las = 1)
mtext('(b)', side = 2, line = 1.5, at = 1.2, las = 1)

#dev.off()



Figure 4

## Plot posterior density

#png("figure4.png", width = 7, height = 7, units = "in", res = 360)
#pdf("figure4.pdf")


par(las = 1)

plot(NA, xlim = c(-1, 1.5), ylim = c(0, 3), axes = FALSE, xlab = '', ylab = '')
grid()

par(new = TRUE)
plot(density(mcmc.model), xlim = c(-1, 1.5))
box()

## add whiskers to show HPD

h <- 0-par()$usr[3]
lines(x = hpd.r, y = rep(-h/2, 2))
lines(x = rep(hpd.r[1], 2), y = c(-0.3, -0.7)*h)
lines(x = rep(hpd.r[2], 2), y = c(-0.3, -0.7)*h)

#dev.off()



Sensitivity analysis - Added by Derek Benson

## Create sensitivity analysis object using 'samp' method

samp <- samp_phylm(Average_Telomere_Length_kb ~ log.lifespan, phy = full_data_tree, data = pruned_dat, n.sim = 1000, track = FALSE)

summary(samp)
>   % Species Removed % Significant Intercepts Mean Intercept Change (%) Mean sDIFintercept % Significant Estimates
> 1                10                     99.9                    3.3375        -0.09634809                    99.8
> 2                20                     99.4                    5.4386        -0.20865625                    98.9
> 3                30                     99.1                    7.2334        -0.32696877                    97.3
> 4                40                     97.9                    9.3075        -0.42784789                    91.1
> 5                50                     98.3                   11.4862        -0.49835594                    81.3
>   Mean Estimate Change (%) Mean sDIFestimate
> 1                   6.0234       0.011289138
> 2                   9.9227       0.027871516
> 3                  12.9821       0.052113105
> 4                  16.9222       0.099651154
> 5                  20.9871       0.006694003
sensi_plot(samp, graphs = 1, param = "estimate")

sensi_plot(samp, graphs = 2, param = "estimate")

sensi_plot(samp, graphs = 3, param = "estimate")

sensi_plot(samp, graphs = 4, param = "estimate")